Series-apple

Acciones de APPLE

Cierre de las acciones de APPLE desde el 2015-05-27 al 2020-05-10 medidas diariamente.

AAPL <- read_csv("AAPL.csv")

## New names:
## Rows: 1258 Columns: 15
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," chr
## (1): symbol dbl (13): ...1, close, high, low, open, volume, adjClose, adjHigh,
## adjLow, ... dttm (1): date
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## • `` -> `...1`

apple2 <- xts(x = AAPL[,"close"],
                  order.by = AAPL$date)

TSstudio::ts_plot(apple2,title="",slider=TRUE)

Del grafico podemos observar que se tiene tendencia que podría ser lineal además de que su varianza va aumentando a traves del tiempo, no parece tener estacionalidad.

# REGULARIDAD DE LA SERIE
is.regular(apple2,strict = TRUE)

## [1] FALSE

is.regular(apple2,strict = FALSE)

## [1] TRUE

Decimos entonces que la serie no es estrictamente regular puesto que no obtiene valores los fines de semana.

Estabilización de la varianza

MASS::boxcox(lm(apple2 ~ 1),seq(-2, 3, length = 50))  

forecast::BoxCox.lambda(apple2, method ="loglik", lower = -1, upper = 3)

## [1] -0.35

El valor de \(\lambda = -0.35\) en este caso, usaremos a box-cox con \(\lambda^{-1}(\mu_t^\lambda -1)\).

plot(BoxCox(apple2,lambda = -0.35))

lambda <- -0.35
B_apple<- lambda^{-1}*((apple2^lambda)-1)
plot(B_apple)

# Objeto ts con estabilización de la varianza
B_apple_ts <- as.ts(coredata(B_apple))

Se ve un cambio en la pequeño en la varianza, tras la transformación.

Estimación y eliminación de tendencia

Tomando los datos con la varianza estabilizada, ajustamos una regresión lineal a la serie de tiempo.

#Serie de tiempo en objeto ts
B_apple_ts <- as.ts(coredata(B_apple))

summary(fit_apple <- lm(B_apple_ts~time(B_apple_ts), na.action=NULL))

## 
## Call:
## lm(formula = B_apple_ts ~ time(B_apple_ts), na.action = NULL)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.054988 -0.017148  0.003503  0.014529  0.058230 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      2.282e+00  1.221e-03 1868.68   <2e-16 ***
## time(B_apple_ts) 1.393e-04  1.680e-06   82.89   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.02164 on 1256 degrees of freedom
## Multiple R-squared:  0.8455, Adjusted R-squared:  0.8453 
## F-statistic:  6871 on 1 and 1256 DF,  p-value: < 2.2e-16

plot(B_apple_ts, ylab="Número de Pasajeros en escala logarítmica") 
abline(fit_apple,col = "red") # Se añade la recta ajusta

Se concluye que la acción tiene una tendincia lineal creciente.

Por otro lado con por medio de un enfoque no parametrico, por medio de un enfoque de STL podemos notar una tendencia no lineal.

indice_bapple <- as.Date(tk_index(B_apple))
df_bapple <- data.frame(Fecha=indice_bapple, as.matrix(B_apple))
tibble_lapple <- tibble(df_bapple)

duplicates(tibble_lapple, key = NULL, index=Fecha)   

## # A tibble: 0 × 2
## # ℹ 2 variables: Fecha <date>, close <dbl>

###Ajuste STL moviendo los parámetros
tibble_lapple%>%mutate(Lapple_ajus=smooth_vec(close,span = 0.75, degree = 2))

## # A tibble: 1,258 × 3
##    Fecha      close Lapple_ajus
##    <date>     <dbl>       <dbl>
##  1 2015-05-27  2.34        2.32
##  2 2015-05-28  2.34        2.32
##  3 2015-05-29  2.34        2.32
##  4 2015-06-01  2.34        2.32
##  5 2015-06-02  2.34        2.32
##  6 2015-06-03  2.34        2.32
##  7 2015-06-04  2.34        2.32
##  8 2015-06-05  2.34        2.32
##  9 2015-06-08  2.33        2.32
## 10 2015-06-09  2.33        2.32
## # ℹ 1,248 more rows

tibble_lapple%>%mutate(Lapple_ajus=smooth_vec(close,span = 0.75, degree = 2))%>%
  ggplot(aes(Fecha, close)) +
  geom_line() +
  geom_line(aes(y = Lapple_ajus), color = "red")

Por otro lado la descomposición por medio de STL, tenemos:

tsibble_lapple1<-as_tsibble(B_apple_ts)
str(tsibble_lapple1)

## tbl_ts [1,258 × 2] (S3: tbl_ts/tbl_df/tbl/data.frame)
##  $ index: num [1:1258] 1 2 3 4 5 6 7 8 9 10 ...
##  $ value: num [1:1258] 2.34 2.34 2.34 2.34 2.34 ...
##  - attr(*, "key")= tibble [1 × 1] (S3: tbl_df/tbl/data.frame)
##   ..$ .rows: list<int> [1:1] 
##   .. ..$ : int [1:1258] 1 2 3 4 5 6 7 8 9 10 ...
##   .. ..@ ptype: int(0) 
##  - attr(*, "index")= chr "index"
##   ..- attr(*, "ordered")= logi TRUE
##  - attr(*, "index2")= chr "index"
##  - attr(*, "interval")= interval [1:1] 1
##   ..@ .regular: logi TRUE

tsibble_lapple1 %>%
  model(
    STL(value ~ trend() +
          season(window = "periodic"),
        robust = TRUE)) %>%
  components() %>%
  autoplot()

Para eliminar la tendencia vamos a utilizar diferencia ordinaria:

###Diferenciando basado en el objeto tibble
tibble_lapple%>%mutate(diff_Lapple=close-lag(close))%>%plot_time_series(Fecha,diff_Lapple)

tibble_lapple<-tibble_lapple%>%mutate(diff_Lapple=close-lag(close))

## con ts
dlapple<-diff(B_apple_ts)

Relaciones no lineales

En un inicio buscaremos relaciones de la variable respecto al tiempo por sus retargos:

par(mar = c(3,2,3,2))
astsa::lag1.plot(
  dlapple,15,corr=T)

Con 16 retargos no hay relaciones significativas observadas en el tiempo.

acf(dlapple, 48, main="Serie diferenciada y con BoxCox de Apple",na.action = na.pass)

pacf(dlapple, 48, na.action = na.pass)

En cuanto a el AMI, tenemos que las variables no explican mucho una de la otra dado el tiempo.

nonlinearTseries::mutualInformation(dlapple,lag.max = 100,n.partitions = 50,units = "Bits",do.plot = TRUE) #c("Nats", "Bits", "Bans")

## $time.lag
##   [1]   0   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17
##  [19]  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35
##  [37]  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53
##  [55]  54  55  56  57  58  59  60  61  62  63  64  65  66  67  68  69  70  71
##  [73]  72  73  74  75  76  77  78  79  80  81  82  83  84  85  86  87  88  89
##  [91]  90  91  92  93  94  95  96  97  98  99 100
## 
## $mutual.information
##   [1] 3.8934615 0.4295910 0.4576206 0.4372052 0.4499625 0.4252671 0.4401849
##   [8] 0.4308357 0.4556016 0.4211155 0.4025746 0.4559729 0.4056106 0.4349552
##  [15] 0.4152535 0.3921074 0.4047687 0.4124704 0.4148666 0.3994185 0.4236338
##  [22] 0.3844587 0.3913508 0.4113851 0.3858057 0.3848515 0.3927395 0.4100664
##  [29] 0.3578133 0.3867007 0.4013168 0.4061538 0.3900069 0.3759358 0.3975308
##  [36] 0.3803090 0.3675074 0.3664258 0.3742754 0.3487143 0.3930242 0.3657913
##  [43] 0.3582033 0.3525726 0.3350029 0.3386703 0.3299502 0.3089034 0.3473679
##  [50] 0.3238666 0.3167121 0.3055256 0.2978825 0.3163190 0.2774472 0.3291654
##  [57] 0.2756802 0.2751606 0.2777094 0.2700078 0.2877226 0.2996385 0.2595825
##  [64] 0.2677796 0.2715066 0.2514312 0.2681325 0.2825981 0.2379677 0.2725874
##  [71] 0.2930380 0.3142999 0.2732301 0.2774603 0.2877546 0.2808847 0.2900511
##  [78] 0.2815907 0.2784265 0.2941168 0.2729975 0.2723005 0.3056265 0.2840801
##  [85] 0.2944494 0.2851070 0.3091139 0.2859998 0.2968052 0.3220082 0.2753343
##  [92] 0.2864599 0.3011603 0.3108419 0.3090950 0.3150884 0.3148816 0.3118721
##  [99] 0.2983499 0.3159033 0.3080346
## 
## $units
## [1] "Bits"
## 
## $n.partitions
## [1] 50
## 
## attr(,"class")
## [1] "mutualInf"

Deteccion de estacionalidad

Realizamos un acercamiento a la detección de la estacionalidad.

El grafico muestra el comportamiento medio de las acciones diariamente en los años.

tibble_lapple %>%na.omit()|>
  mutate(
    DIA = str_c("", as.character(lubridate::day(Fecha)))
  ) %>%
  plot_time_series(
    .date_var = Fecha,
    .value = diff_Lapple,
    .facet_vars = DIA,
    .facet_ncol = 4, 
    .color_var = DIA, 
    .facet_scales = "fixed",
    .interactive = FALSE,
    .legend_show = FALSE,
    .smooth = FALSE
  )

Donde se observa que dados los diferentes días el comportamiento medio no varia mucho uno del otro.

Se realiza un Boxplot buscando posible estacionalidad dado el dia de la semana, la semana y el mes.

## Al hacerla por dias de la semana no hay muchas diferencias de medias
tibble_lapple%>%na.omit()%>%
  plot_seasonal_diagnostics(.date_var = Fecha,.value = diff_Lapple,.feature_set = c("wday.lbl"),.geom="boxplot")

## Al hacerla por semana no hay diferencias de medias pero no son muy marcadas
tibble_lapple%>%na.omit()%>%
  plot_seasonal_diagnostics(.date_var = Fecha,.value = diff_Lapple,.feature_set = c("week"),.geom="boxplot")

## Al hacerla por mes no hay muchas diferencias de medias
tibble_lapple%>%na.omit()%>%
  plot_seasonal_diagnostics(.date_var = Fecha,.value = diff_Lapple,.feature_set = c("month.lbl"),.geom="boxplot")

De los anteriores graficos se observa que no hay una diferencia significativa de de medias en ninguno de los tres casos mencionados.

Por ultimo al realizar graficos de densidad por kernell llegamos a una conclusión similar dado el día.

## Por kernell las densidades diarias no varian mucho más de las otras
ggplot(tibble_lapple %>%na.omit()|>
         mutate(
           DIA = str_c("", as.character(lubridate::day(Fecha)))
         ), aes(x = diff_Lapple)) +
  geom_density(aes(fill = DIA)) +
  ggtitle("LosPass - Estimación de la densidad vía Kernel por mes") +
  facet_grid(rows = vars(as.factor(DIA)))

Suavizamiento Exponencial

Con lo anterior, se va a trabajar con el objeto de serie de tiempo tras la transformación de Box-Cox.

Los graficos presentan la modelación de la tendencia con paramentros \(\alpha\), el nivel, y un pronostico de 24 observaciones. No se hace uso de \(\gamma\) puesto que no se identifico una componente de estacionalidad en la serie presentada. (SEASONAL es el default)

HWAP_inic=stats::HoltWinters(B_apple_ts,alpha=NULL,beta=FALSE,gamma=FALSE)
plot(HWAP_inic)

plot(forecast::forecast(HWAP_inic,h=24,level =0.95))

De lo anterior podemos evidenciar que el modelo tiene un buen ajuste frente a los datos presentados, en cuanto a su pronostico parece no continua con la tendencia además sus intervalos parecen ser un poco grandes.

HWAP=stats::HoltWinters(B_apple_ts,gamma=FALSE,seasonal="additive")
plot(HWAP)

plot(forecast::forecast(HWAP,h=24,level =0.95))

Con el modelo anterior tenemos el parametro de \(\alpha\) y \(\beta\), la pendiente, utilizando un modelo de tipo aditivo, donde notamos un buen ajuste frente a los datos presentados, frente al pronostico podemos ver que se mantiene la tendencia y que los intervalos de confianza no varian de forma significativa frente al modelo anterior.

ajustados=fitted(HWAP)
plot(ajustados)

HWAP

## Holt-Winters exponential smoothing with trend and without seasonal component.
## 
## Call:
## stats::HoltWinters(x = B_apple_ts, gamma = FALSE, seasonal = "additive")
## 
## Smoothing parameters:
##  alpha: 0.9221333
##  beta : 0.004817684
##  gamma: FALSE
## 
## Coefficients:
##           [,1]
## a 2.4772099306
## b 0.0002218824

Haciendo uso de un tsibble

ajustepass=tsibble_lapple1%>%
  model(ETS(value~ error("A")+trend("A")+season("N")))

pronostico=ajustepass%>%
  fabletools::forecast(h=12)
pronostico

## # A fable: 12 x 4 [1]
## # Key:     .model [1]
##    .model                                            index           value .mean
##    <chr>                                             <dbl>          <dist> <dbl>
##  1 "ETS(value ~ error(\"A\") + trend(\"A\") + seaso…  1259   N(2.5, 9e-06)  2.48
##  2 "ETS(value ~ error(\"A\") + trend(\"A\") + seaso…  1260 N(2.5, 1.7e-05)  2.48
##  3 "ETS(value ~ error(\"A\") + trend(\"A\") + seaso…  1261 N(2.5, 2.4e-05)  2.48
##  4 "ETS(value ~ error(\"A\") + trend(\"A\") + seaso…  1262 N(2.5, 3.2e-05)  2.48
##  5 "ETS(value ~ error(\"A\") + trend(\"A\") + seaso…  1263   N(2.5, 4e-05)  2.48
##  6 "ETS(value ~ error(\"A\") + trend(\"A\") + seaso…  1264 N(2.5, 4.7e-05)  2.48
##  7 "ETS(value ~ error(\"A\") + trend(\"A\") + seaso…  1265 N(2.5, 5.5e-05)  2.48
##  8 "ETS(value ~ error(\"A\") + trend(\"A\") + seaso…  1266 N(2.5, 6.3e-05)  2.48
##  9 "ETS(value ~ error(\"A\") + trend(\"A\") + seaso…  1267   N(2.5, 7e-05)  2.48
## 10 "ETS(value ~ error(\"A\") + trend(\"A\") + seaso…  1268 N(2.5, 7.8e-05)  2.48
## 11 "ETS(value ~ error(\"A\") + trend(\"A\") + seaso…  1269 N(2.5, 8.6e-05)  2.48
## 12 "ETS(value ~ error(\"A\") + trend(\"A\") + seaso…  1270 N(2.5, 9.3e-05)  2.48

pronostico%>%autoplot(tsibble_lapple1)+geom_line(aes(y=.fitted),col="#D55E00",data=augment(ajustepass))+labs(y=" ",title="Pronóstico u ajustados")+guides(colour="none")

POR ULTIMO

modelos<-tsibble_lapple1%>%
   model(ets=ETS(value~ error("A")+trend("A")+season("N")),
         stl=decomposition_model(STL(value ~ trend(window = 13) +
                   season(window = "periodic"),
    robust = TRUE),NAIVE(season_adjust)))
 modelos 

## # A mable: 1 x 2
##            ets                       stl
##        <model>                   <model>
## 1 <ETS(A,A,N)> <STL decomposition model>

 modelos%>%fabletools::forecast(h=12)%>%
   autoplot(tsibble_lapple1)

Arboles de desición

Se hace uso de la serie sin tndencia.

data = pd.read_csv('tibble_lapple.txt', skiprows=1, delimiter=',', header = 0, usecols =[1,3],
                   names =["Fecha", "diff lapple"], dtype={"diff lapple": np.float64})
dlapple=data.set_index("Fecha")
t_dlapple_v = pd.DataFrame(data['diff lapple'].values,index=data['Fecha'])
t_dlapple=dlapple['diff lapple']
plt.plot(t_dlapple)
plt.title('Cierre acciones de apple sin tendencia') 

Tenemos la función de autocorrelación y aoutocorrelación parcial

from statsmodels.graphics.tsaplots import plot_acf
plot_acf(t_dlapple,lags=100)
pyplot.show()

from statsmodels.graphics.tsaplots import plot_pacf
plot_pacf(t_dlapple,lags=100,method='ywm',alpha=0.01)
pyplot.show()

Además de esto, utilizaremos el periodograma para determinar las covariables que se van a utilizar

f, Pxx_den=sp.signal.periodogram(t_dlapple.values)
plt.plot(f, Pxx_den)
max_index_value = np.argmax(Pxx_den, axis=0)
print(max_index_value)

## 12

frecuencia_max=f[max_index_value]
print(frecuencia_max)

## 0.00954653937947494

1/frecuencia_max

## 104.75

#print(heapq.nlargest(10, range(len(Pxx_den)), key=Pxx_den.__getitem__))

#print(1/f[4])
#print(1/f[5])
#print(1/f[9])
#print(1/f[7])
#print(1/f[15])
#print(1/f[11])

f_welch, Pxx_den_welch=sp.signal.welch(t_dlapple.values)
plt.plot(f_welch, Pxx_den_welch)

max_index_value_welch = np.argmax(Pxx_den_welch, axis=0)
print(max_index_value_welch)

## 109

frecuencia_max_welch=f_welch[max_index_value_welch]
print(frecuencia_max_welch)

## 0.42578125

1/frecuencia_max_welch

## 2.3486238532110093

import heapq
print(heapq.nlargest(10, range(len(Pxx_den_welch)), key=Pxx_den_welch.__getitem__))

## [109, 119, 40, 108, 30, 110, 2, 73, 31, 92]

print(1/f_welch[1])

## 256.0

print(1/f_welch[2])

## 128.0

print(1/f_welch[3])

## 85.33333333333333

Dado lo anterior, utilizaremos como covariables de los primeros 12 retardos cada dos y los retardos del 120 al 90

var_rez = pd.DataFrame()

for i in range(12,0,-2):
    var_rez[['t-'+str(i)]] = t_dlapple_v.shift(i)

for i in range(190,120,-1):
    var_rez[['t-'+str(i)]] = t_dlapple_v.shift(i)

#for i in range(560,530,-2):
 #   var_rez[['t-'+str(i)]] = t_dlapple_v.shift(i)
 
var_rez['t'] = t_dlapple_v.values
var_rezl = var_rez[190:]
print(var_rezl)

##                 t-12      t-10       t-8  ...     t-122     t-121         t
## Fecha                                     ...                              
## 2016-02-29 -0.001548  0.000630  0.003061  ...  0.008104 -0.003400 -0.000459
## 2016-03-01 -0.001237  0.005642 -0.003857  ... -0.003400 -0.001934  0.007810
## 2016-03-02  0.000630  0.003061 -0.000463  ... -0.001934  0.005283  0.000435
## 2016-03-03  0.005642 -0.003857  0.001760  ...  0.005283 -0.003733  0.001474
## 2016-03-04  0.003061 -0.000463 -0.004632  ... -0.003733  0.004176  0.002924
## ...              ...       ...       ...  ...       ...       ...       ...
## 2020-05-18  0.002867  0.001929  0.001396  ... -0.001662 -0.000639  0.003123
## 2020-05-19 -0.002227  0.002035  0.001394  ... -0.000639 -0.000125 -0.000775
## 2020-05-20  0.001929  0.001396  0.002806  ... -0.000125  0.002469  0.002569
## 2020-05-21  0.002035  0.001394  0.002091  ...  0.002469 -0.001112 -0.000996
## 2020-05-22  0.001396  0.002806 -0.001538  ... -0.001112  0.001890  0.000854
## 
## [1067 rows x 77 columns]

Se realiza la división de datos para entrenamiento, validación y prueba para la “respuesta” y para sus covariables

dlapple_split = var_rezl.values
X1= dlapple_split[:, 0:-1]  
y1 =dlapple_split[:,-1]  

Y1 = y1
print('Complete Observations for Target after Supervised configuration: %d' %len(Y1))

## Complete Observations for Target after Supervised configuration: 1067

traintarget_size = int(len(Y1) * 0.70) 
valtarget_size = int(len(Y1) * 0.10)# Set split
testtarget_size = int(len(Y1) * 0.20)# Set split
print(traintarget_size,valtarget_size,testtarget_size)

## 746 106 213

print('Train + Validation + Test: %d' %(traintarget_size+valtarget_size+testtarget_size))

## Train + Validation + Test: 1065

Y1 = y1
traintarget_size = int(len(Y1) * 0.70) 
valtarget_size = int(len(Y1) * 0.10)+1# Set split
testtarget_size = int(len(Y1) * 0.20)+1# Set split
train_target, val_target,test_target = Y1[0:traintarget_size],Y1[(traintarget_size):(traintarget_size+valtarget_size)] ,Y1[(traintarget_size+valtarget_size):len(Y1)]

print('Observations for Target: %d' % (len(Y1)))

## Observations for Target: 1067

print('Training Observations for Target: %d' % (len(train_target)))

## Training Observations for Target: 746

print('Validation Observations for Target: %d' % (len(val_target)))

## Validation Observations for Target: 107

print('Test Observations for Target: %d' % (len(test_target)))

## Test Observations for Target: 214

trainfeature_size = int(len(X1) * 0.70)
valfeature_size = int(len(X1) * 0.10)+1# Set split
testfeature_size = int(len(X1) * 0.20)# Set split
train_feature, val_feature,test_feature = X1[0:traintarget_size],X1[(traintarget_size):(traintarget_size+valtarget_size)] ,X1[(traintarget_size+valtarget_size):len(Y1)]

print('Observations for Feature: %d' % (len(X1)))

## Observations for Feature: 1067

print('Training Observations for Feature: %d' % (len(train_feature)))

## Training Observations for Feature: 746

print('Validation Observations for Feature: %d' % (len(val_feature)))

## Validation Observations for Feature: 107

print('Test Observations for Feature: %d' % (len(test_feature)))

## Test Observations for Feature: 214

Se aplica un arbol de regresión sin profundidad maxima estipulada

from sklearn.tree import DecisionTreeRegressor

decision_tree_PM25 = DecisionTreeRegressor()  # max-depth not set
decision_tree_PM25.fit(train_feature, train_target)

DecisionTreeRegressor()

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

print("Coeficiente R2 sobre el conjunto de entrenamiento:",decision_tree_PM25.score(train_feature, train_target))

## Coeficiente R2 sobre el conjunto de entrenamiento: 1.0

print("Coeficiente R2 sobre el conjunto de Validación:",decision_tree_PM25.score(val_feature,val_target))  # predictions are horrible if negative value, no relationship if 0

## Coeficiente R2 sobre el conjunto de Validación: -1.7175114296718554

print("el RECM sobre validación es:",(((decision_tree_PM25.predict(val_feature)-val_target)**2).mean()) )

## el RECM sobre validación es: 1.4301750781635833e-05

Note que si bien se tiene un ajuste perfecto en entrenamiento, en validación el \(R^2\) tiene un valor negativo, lo que implica un sobreajuste sobre los datos de entrenamiento.

for d in [2, 3, 4, 5,6,8,9,10,11,12]:
    decision_tree_PM25 = DecisionTreeRegressor(max_depth=d)
    decision_tree_PM25.fit(train_feature, train_target)

    print('max_depth=', str(d))
    print("Coeficiente R2 sobre el conjunto de entrenamiento:",decision_tree_PM25.score(train_feature, train_target))
    print("Coeficiente R2 sobre el conjunto de validación:",decision_tree_PM25.score(val_feature, val_target), '\n')  # You want the test score to be positive and high
    print("el RECM sobre el conjunto de validación es:",sklearn.metrics.mean_squared_error(decision_tree_PM25.predict(val_feature),val_target, squared=False))

DecisionTreeRegressor(max_depth=12)

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

Se usa una grilla de valores para evaluar un posible hiperparametro para la profundidad, en este caso esta será de 2 usando el coeficiente de \(R^2\). Ya definido esto, se reentrena el modelo con los datos de entrenamiento y validación.

train_val_feature=np.concatenate((train_feature,val_feature),axis=0)
train_val_target=np.concatenate((train_target,val_target),axis=0)

from matplotlib import pyplot as plt

decision_tree_PM25 = DecisionTreeRegressor(max_depth=2)  # fill in best max depth here
decision_tree_PM25.fit(train_val_feature, train_val_target)

DecisionTreeRegressor(max_depth=2)

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

train_val_prediction = decision_tree_PM25.predict(train_val_feature)
test_prediction = decision_tree_PM25.predict(test_feature)

plt.scatter(train_val_prediction, train_val_target, label='train')  # blue
plt.scatter(test_prediction, test_target, label='test')  # orange
plt.show()

print("Raíz de la Pérdida cuadrática Entrenamiento:",sklearn.metrics.mean_squared_error( train_val_prediction, train_val_target,squared=False))

## Raíz de la Pérdida cuadrática Entrenamiento: 0.002533460453410802

print("Raíz de la Pérdida cuadrática Prueba:",sklearn.metrics.mean_squared_error(test_prediction, test_target,squared=False))

## Raíz de la Pérdida cuadrática Prueba: 0.0038384586318421744

De lo anterior podemos concluir que lo modelado tiente a tomar en entrenamiento datos en la media comparado con su valor real.

A continuación se presenta un gráfico de las predicciones del entrenamiento y de prueba, donde ambas tienden a la media.

from matplotlib import pyplot as plt

decision_tree_PM25_prun_mincost = DecisionTreeRegressor(max_depth=2)  # fill in best max depth here
decision_tree_PM25_prun_mincost.fit(train_val_feature, train_val_target)

DecisionTreeRegressor(max_depth=2)

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

train_val_prediction_prun_mincost = decision_tree_PM25.predict(train_val_feature)
test_prediction_prun_mincost = decision_tree_PM25.predict(test_feature)

plt.scatter(train_val_prediction_prun_mincost, train_val_target, label='train')  # blue
plt.scatter(test_prediction_prun_mincost, test_target, label='test')  # orange
plt.show()

indicetrian_val_test=var_rezl.index
print(indicetrian_val_test.size)  ###Tamaño del índice

## 1067

indicetrain_val=indicetrian_val_test[0:214]
indicetest=indicetrian_val_test[214:853]

targetjoint=np.concatenate((train_val_target,test_target))

predictionjoint=np.concatenate((train_val_prediction,test_prediction))

d = {'observado': targetjoint, 'Predicción': predictionjoint}
ObsvsPred=pd.DataFrame(data=d,index=indicetrian_val_test)

En la siguiente gráfica se muestra el valor real frete a su predicción.

plt.plot(ObsvsPred)